home *** CD-ROM | disk | FTP | other *** search
/ Info-Mac 11 / Info-Mac_XI_Disc_1.cdr_ / Info-Mac XI Disc 1.cdr / Programs / Science & Math / MacEspresso 1.0 / espresso / set.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-12-03  |  17.2 KB  |  811 lines  |  [TEXT/R*ch]
  1. /*
  2.  *   set.c -- routines for maniuplating sets and set families
  3.  */
  4.  
  5. /* LINTLIBRARY */
  6.  
  7. #include "espresso.h"
  8. static pset_family set_family_garbage = NULL;
  9.  
  10. static int intcpy(d, s, n)
  11. register unsigned int *d, *s;
  12. register long n;
  13. {
  14.     register int i;
  15.     for(i = 0; i < n; i++) {
  16.     *d++ = *s++;
  17.     }
  18. }
  19.  
  20.  
  21. /* bit_index -- find first bit (from LSB) in a word (MSB=bit n, LSB=bit 0) */
  22. int bit_index(a)
  23. register unsigned int a;
  24. {
  25.     register int i;
  26.     if (a == 0)
  27.     return -1;
  28.     for(i = 0; (a & 1) == 0; a >>= 1, i++)
  29.     ;
  30.     return i;
  31. }
  32.  
  33.  
  34. /* set_ord -- count number of elements in a set */
  35. int set_ord(a)
  36. register pset a;
  37. {
  38.     register int i, sum = 0;
  39.     register unsigned int val;
  40.     for(i = LOOP(a); i > 0; i--)
  41.     if ((val = a[i]) != 0)
  42.         sum += count_ones(val);
  43.     return sum;
  44. }
  45.  
  46. /* set_dist -- distance between two sets (# elements in common) */
  47. int set_dist(a, b)
  48. register pset a, b;
  49. {
  50.     register int i, sum = 0;
  51.     register unsigned int val;
  52.     for(i = LOOP(a); i > 0; i--)
  53.     if ((val = a[i] & b[i]) != 0)
  54.         sum += count_ones(val);
  55.     return sum;
  56. }
  57.  
  58. /* set_clear -- make "r" the empty set of "size" elements */
  59. pset set_clear(r, size)
  60. register pset r;
  61. int size;
  62. {
  63.     register int i = LOOPINIT(size);
  64.     *r = i; do r[i] = 0; while (--i > 0);
  65.     return r;
  66. }
  67.  
  68. /* set_fill -- make "r" the universal set of "size" elements */
  69. pset set_fill(r, size)
  70. register pset r;
  71. register int size;
  72. {
  73.     register int i = LOOPINIT(size);
  74.     *r = i;
  75.     r[i] = ~ (unsigned) 0;
  76.     r[i] >>= i * BPI - size;
  77.     while (--i > 0)
  78.     r[i] = ~ (unsigned) 0;
  79.     return r;
  80. }
  81.  
  82. /* set_copy -- copy set a into set r */
  83. pset set_copy(r, a)
  84. register pset r, a;
  85. {
  86.     register int i = LOOPCOPY(a);
  87.     do r[i] = a[i]; while (--i >= 0);
  88.     return r;
  89. }
  90.  
  91. /* set_and -- compute intersection of sets "a" and "b" */
  92. pset set_and(r, a, b)
  93. register pset r, a, b;
  94. {
  95.     register int i = LOOP(a);
  96.     PUTLOOP(r,i); do r[i] = a[i] & b[i]; while (--i > 0);
  97.     return r;
  98. }
  99.  
  100. /* set_or -- compute union of sets "a" and "b" */
  101. pset set_or(r, a, b)
  102. register pset r, a, b;
  103. {
  104.     register int i = LOOP(a);
  105.     PUTLOOP(r,i); do r[i] = a[i] | b[i]; while (--i > 0);
  106.     return r;
  107. }
  108.  
  109. /* set_diff -- compute difference of sets "a" and "b" */
  110. pset set_diff(r, a, b)
  111. register pset r, a, b;
  112. {
  113.     register int i = LOOP(a);
  114.     PUTLOOP(r,i); do r[i] = a[i] & ~b[i]; while (--i > 0);
  115.     return r;
  116. }
  117.  
  118. /* set_xor -- compute exclusive-or of sets "a" and "b" */
  119. pset set_xor(r, a, b)
  120. register pset r, a, b;
  121. {
  122.     register int i = LOOP(a);
  123. #ifdef IBM_WATC
  124.     PUTLOOP(r,i); do r[i] = (a[i]&~b[i]) | (~a[i]&b[i]); while (--i > 0);
  125. #else
  126.     PUTLOOP(r,i); do r[i] = a[i] ^ b[i]; while (--i > 0);
  127. #endif
  128.     return r;
  129. }
  130.  
  131. /* set_merge -- compute "a" & "mask" | "b" & ~ "mask" */
  132. pset set_merge(r, a, b, mask)
  133. register pset r, a, b, mask;
  134. {
  135.     register int i = LOOP(a);
  136.     PUTLOOP(r,i); do r[i] = (a[i]&mask[i]) | (b[i]&~mask[i]); while (--i > 0);
  137.     return r;
  138. }
  139.  
  140. /* set_andp -- compute intersection of sets "a" and "b" , TRUE if nonempty */
  141. bool set_andp(r, a, b)
  142. register pset r, a, b;
  143. {
  144.     register int i = LOOP(a);
  145.     register unsigned int x = 0;
  146.     PUTLOOP(r,i); do {r[i] = a[i] & b[i]; x |= r[i];} while (--i > 0);
  147.     return x != 0;
  148. }
  149.  
  150. /* set_orp -- compute union of sets "a" and "b" , TRUE if nonempty */
  151. bool set_orp(r, a, b)
  152. register pset r, a, b;
  153. {
  154.     register int i = LOOP(a);
  155.     register unsigned int x = 0;
  156.     PUTLOOP(r,i); do {r[i] = a[i] | b[i]; x |= r[i];} while (--i > 0);
  157.     return x != 0;
  158. }
  159.  
  160. /* setp_empty -- check if the set "a" is empty */
  161. bool setp_empty(a)
  162. register pset a;
  163. {
  164.     register int i = LOOP(a);
  165.     do if (a[i]) return FALSE; while (--i > 0);
  166.     return TRUE;
  167. }
  168.  
  169. /* setp_full -- check if the set "a" is the full set of "size" elements */
  170. bool setp_full(a, size)
  171. register pset a;
  172. register int size;
  173. {
  174.     register int i = LOOP(a);
  175.     register unsigned int test;
  176.     test = ~ (unsigned) 0;
  177.     test >>= i * BPI - size;
  178.     if (a[i] != test)
  179.     return FALSE;
  180.     while (--i > 0)
  181.     if (a[i] != (~(unsigned) 0))
  182.         return FALSE;
  183.     return TRUE;
  184. }
  185.  
  186. /* setp_equal -- check if the set "a" equals set "b" */
  187. bool setp_equal(a, b)
  188. register pset a, b;
  189. {
  190.     register int i = LOOP(a);
  191.     do if (a[i] != b[i]) return FALSE; while (--i > 0);
  192.     return TRUE;
  193. }
  194.  
  195. /* setp_disjoint -- check if intersection of "a" and "b" is empty */
  196. bool setp_disjoint(a, b)
  197. register pset a, b;
  198. {
  199.     register int i = LOOP(a);
  200.     do if (a[i] & b[i]) return FALSE; while (--i > 0);
  201.     return TRUE;
  202. }
  203.  
  204. /* setp_implies -- check if "a" implies "b" ("b" contains "a") */
  205. bool setp_implies(register pset a,register pset b)
  206. {
  207.     register int i = LOOP(a);
  208.     do if (a[i] & ~b[i]) return FALSE; while (--i > 0);
  209.     return TRUE;
  210. }
  211.  
  212. /* sf_or -- form the "or" of all sets in a set family */
  213. pset sf_or(A)
  214. pset_family A;
  215. {
  216.     register pset or, last, p;
  217.  
  218.     or = set_new(A->sf_size);
  219.     foreach_set(A, last, p)
  220.     INLINEset_or(or, or, p);
  221.     return or;
  222. }
  223.  
  224. /* sf_and -- form the "and" of all sets in a set family */
  225. pset sf_and(A)
  226. pset_family A;
  227. {
  228.     register pset and, last, p;
  229.  
  230.     and = set_fill(set_new(A->sf_size), A->sf_size);
  231.     foreach_set(A, last, p)
  232.     INLINEset_and(and, and, p);
  233.     return and;
  234. }
  235.  
  236. /* sf_active -- make all members of the set family active */
  237. pset_family sf_active(A)
  238. pset_family A;
  239. {
  240.     register pset p, last;
  241.     foreach_set(A, last, p) {
  242.     SET(p, ACTIVE);
  243.     }
  244.     A->active_count = A->count;
  245.     return A;
  246. }
  247.  
  248.  
  249. /* sf_inactive -- remove all inactive cubes in a set family */
  250. pset_family sf_inactive(A)
  251. pset_family A;
  252. {
  253.     register pset p, last, pdest;
  254.  
  255.     pdest = A->data;
  256.     foreach_set(A, last, p) {
  257.     if (TESTP(p, ACTIVE)) {
  258.         if (pdest != p) {
  259.         INLINEset_copy(pdest, p);
  260.         }
  261.         pdest += A->wsize;
  262.     } else {
  263.         A->count--;
  264.     }
  265.     }
  266.     return A;
  267. }
  268.  
  269.  
  270. /* sf_copy -- copy a set family */
  271. pset_family sf_copy(R, A)
  272. pset_family R, A;
  273. {
  274.     R->sf_size = A->sf_size;
  275.     R->wsize = A->wsize;
  276. /*R->capacity = A->count;*/
  277. /*R->data = REALLOC(unsigned int, R->data, (long) R->capacity * R->wsize);*/
  278.     R->count = A->count;
  279.     R->active_count = A->active_count;
  280.     intcpy(R->data, A->data, (long) A->wsize * A->count);
  281.     return R;
  282. }
  283.  
  284.  
  285. /* sf_join -- join A and B into a single set_family */
  286. pset_family sf_join(A, B)
  287. pset_family A, B;
  288. {
  289.     pset_family R;
  290.     long asize = A->count * A->wsize;
  291.     long bsize = B->count * B->wsize;
  292.  
  293.     if (A->sf_size != B->sf_size) fatal("sf_join: sf_size mismatch");
  294.     R = sf_new(A->count + B->count, A->sf_size);
  295.     R->count = A->count + B->count;
  296.     R->active_count = A->active_count + B->active_count;
  297.     intcpy(R->data, A->data, asize);
  298.     intcpy(R->data + asize, B->data, bsize);
  299.     return R;
  300. }
  301.  
  302.  
  303. /* sf_append -- append the sets of B to the end of A, and dispose of B */
  304. pset_family sf_append(A, B)
  305. pset_family A, B;
  306. {
  307.     long asize = A->count * A->wsize;
  308.     long bsize = B->count * B->wsize;
  309.  
  310.     if (A->sf_size != B->sf_size) fatal("sf_append: sf_size mismatch");
  311.     A->capacity = A->count + B->count;
  312.     A->data = REALLOC(unsigned int, A->data, (long) A->capacity * A->wsize);
  313.     intcpy(A->data + asize, B->data, bsize);
  314.     A->count += B->count;
  315.     A->active_count += B->active_count;
  316.     sf_free(B);
  317.     return A;
  318. }
  319.  
  320.  
  321. /* sf_new -- allocate "num" sets of "size" elements each */
  322. pset_family sf_new(num, size)
  323. int num, size;
  324. {
  325.     pset_family A;
  326.     if (set_family_garbage == NULL) {
  327.     A = ALLOC(set_family_t, 1);
  328.     } else {
  329.     A = set_family_garbage;
  330.     set_family_garbage = A->next;
  331.     }
  332.     A->sf_size = size;
  333.     A->wsize = SET_SIZE(size);
  334.     A->capacity = num;
  335.     A->data = ALLOC(unsigned int, (long) A->capacity * A->wsize);
  336.     A->count = 0;
  337.     A->active_count = 0;
  338.     return A;
  339. }
  340.  
  341.  
  342. /* sf_save -- create a duplicate copy of a set family */
  343. pset_family sf_save(A)
  344. register pset_family A;
  345. {
  346.     return sf_copy(sf_new(A->count, A->sf_size), A);
  347. }
  348.  
  349.  
  350. /* sf_free -- free the storage allocated for a set family */
  351. void sf_free(A)
  352. pset_family A;
  353. {
  354.     FREE(A->data);
  355.     A->next = set_family_garbage;
  356.     set_family_garbage = A;
  357. }
  358.  
  359.  
  360. /* sf_cleanup -- free all of the set families from the garbage list */
  361. void sf_cleanup()
  362. {
  363.     register pset_family p, pnext;
  364.     for(p = set_family_garbage; p != (pset_family) NULL; p = pnext) {
  365.     pnext = p->next;
  366.     FREE(p);
  367.     }
  368.     set_family_garbage = (pset_family) NULL;
  369. }
  370.  
  371.  
  372. /* sf_addset -- add a set to the end of a set family */
  373. pset_family sf_addset(A, s)
  374. pset_family A;
  375. pset s;
  376. {
  377.     register pset p;
  378.  
  379.     if (A->count >= A->capacity) {
  380.     A->capacity = A->capacity + A->capacity/2 + 1;
  381.     A->data = REALLOC(unsigned int, A->data, (long) A->capacity * A->wsize);
  382.     }
  383.     p = GETSET(A, A->count++);
  384.     INLINEset_copy(p, s);
  385.     return A;
  386. }
  387.  
  388. /* sf_delset -- delete a set from a set family */
  389. void sf_delset(A, i)
  390. pset_family A;
  391. int i;
  392. {   (void) set_copy(GETSET(A,i), GETSET(A, --A->count));}
  393.  
  394. /* sf_print -- print a set_family as a set (list the element numbers) */
  395. void sf_print(A)
  396. pset_family A;
  397. {
  398.     char *ps1();
  399.     register pset p;
  400.     register int i;
  401.     foreachi_set(A, i, p)
  402.     printf("A[%d] = %s\n", i, ps1(p));
  403. }
  404.  
  405. /* sf_bm_print -- print a set_family as a bit-matrix */
  406. void sf_bm_print(A)
  407. pset_family A;
  408. {
  409.     char *pbv1();
  410.     register pset p;
  411.     register int i;
  412.     foreachi_set(A, i, p)
  413.     printf("[%4d] %s\n", i, pbv1(p, A->sf_size));
  414. }
  415.  
  416.  
  417. /* sf_write -- output a set family in an unintelligable manner */
  418. void sf_write(fp, A)
  419. FILE *fp;
  420. pset_family A;
  421. {
  422.     register pset p, last;
  423.     fprintf(fp, "%d %d\n", A->count, A->sf_size);
  424.     foreach_set(A, last, p)
  425.     set_write(fp, p);
  426.     (void) fflush(fp);
  427. }
  428.  
  429.  
  430. /* sf_read -- read a set family written by sf_write */
  431. pset_family sf_read(fp)
  432. FILE *fp;
  433. {
  434.     int i, j;
  435.     register pset p, last;
  436.     pset_family A;
  437.  
  438.     (void) fscanf(fp, "%d %d\n", &i, &j);
  439.     A = sf_new(i, j);
  440.     A->count = i;
  441.     foreach_set(A, last, p) {
  442.     (void) fscanf(fp, "%x", p);
  443.     for(j = 1; j <= LOOP(p); j++)
  444.         (void) fscanf(fp, "%x", p+j);
  445.     }
  446.     return A;
  447. }
  448.  
  449.  
  450. /* set_write -- output a set in an unintelligable manner */
  451. void set_write(fp, a)
  452. register FILE *fp;
  453. register pset a;
  454. {
  455.     register int n = LOOP(a), j;
  456.  
  457.     for(j = 0; j <= n; j++) {
  458.     fprintf(fp, "%x ", a[j]);
  459.     if ((j+1) % 8 == 0 && j != n)
  460.         fprintf(fp, "\n\t");
  461.     }
  462.     fprintf(fp, "\n");
  463. }
  464.  
  465.  
  466. /* sf_bm_read -- read a set family written by sf_bm_print (almost) */
  467. pset_family sf_bm_read(fp)
  468. FILE *fp;
  469. {
  470.     int i, j, rows, cols;
  471.     register pset pdest;
  472.     pset_family A;
  473.  
  474.     (void) fscanf(fp, "%d %d\n", &rows, &cols);
  475.     A = sf_new(rows, cols);
  476.     for(i = 0; i < rows; i++) {
  477.     pdest = GETSET(A, A->count++);
  478.     (void) set_clear(pdest, A->sf_size);
  479.     for(j = 0; j < cols; j++) {
  480.         switch(getc(fp)) {
  481.         case '0':
  482.             break;
  483.         case '1':
  484.             set_insert(pdest, j);
  485.             break;
  486.         default:
  487.             fatal("Error reading set family");
  488.         }
  489.     }
  490.     if (getc(fp) != '\n') {
  491.         fatal("Error reading set family (at end of line)");
  492.     }
  493.     }
  494.     return A;
  495. }
  496.  
  497.  
  498.  
  499. /* ps1 -- convert a set into a printable string */
  500. #define largest_string 120
  501. static char s1[largest_string];
  502. char *ps1(a)
  503. register pset a;
  504. {
  505.     register int i, num, l, len = 0, n = NELEM(a);
  506.     char temp[20];
  507.     bool first = TRUE;
  508.  
  509.     s1[len++] = '[';
  510.     for(i = 0; i < n; i++)
  511.     if (is_in_set(a,i)) {
  512.         if (! first)
  513.         s1[len++] = ',';
  514.         first = FALSE; num = i;
  515.         /* Generate digits (reverse order) */
  516.         l = 0; do temp[l++] = num % 10 + '0'; while ((num /= 10) > 0);
  517.         /* Copy them back in correct order */
  518.         do s1[len++] = temp[--l]; while (l > 0);
  519.         if (len > largest_string-15) {
  520.         s1[len++] = '.'; s1[len++] = '.'; s1[len++] = '.';
  521.         break;
  522.         }
  523.     }
  524.  
  525.     s1[len++] = ']';
  526.     s1[len++] = '\0';
  527.     return s1;
  528. }
  529.  
  530. /* pbv1 -- print bit-vector */
  531. char *pbv1(s, n)
  532. pset s;
  533. int n;
  534. {
  535.     register int i;
  536.     for(i = 0; i < n; i++)
  537.     s1[i] = is_in_set(s,i) ? '1' : '0';
  538.     s1[n] = '\0';
  539.     return s1;
  540. }
  541.  
  542.  
  543. /* set_adjcnt -- adjust the counts for a set by "weight" */
  544. void
  545. set_adjcnt(a, count, weight)
  546. register pset a;
  547. register int *count, weight;
  548. {
  549.     register int i, base;
  550.     register unsigned int val;
  551.  
  552.     for(i = LOOP(a); i > 0; ) {
  553.     for(val = a[i], base = --i << LOGBPI; val != 0; base++, val >>= 1) {
  554.         if (val & 1) {
  555.         count[base] += weight;
  556.         }
  557.     }
  558.     }
  559. }
  560.  
  561.  
  562.  
  563. /* sf_count -- perform a column sum over a set family */
  564. int *sf_count(A)
  565. pset_family A;
  566. {
  567.     register pset p, last;
  568.     register int i, base, *count;
  569.     register unsigned int val;
  570.  
  571.     count = ALLOC(int, A->sf_size);
  572.     for(i = A->sf_size - 1; i >= 0; i--) {
  573.     count[i] = 0;
  574.     }
  575.  
  576.     foreach_set(A, last, p) {
  577.     for(i = LOOP(p); i > 0; ) {
  578.         for(val = p[i], base = --i << LOGBPI; val != 0; base++, val >>= 1) {
  579.         if (val & 1) {
  580.             count[base]++;
  581.         }
  582.         }
  583.     }
  584.     }
  585.     return count;
  586. }
  587.  
  588.  
  589. /* sf_count_restricted -- perform a column sum over a set family, restricting
  590.  * to only the columns which are in r; also, the columns are weighted by the
  591.  * number of elements which are in each row
  592.  */
  593. int *sf_count_restricted(A, r)
  594. pset_family A;
  595. register pset r;
  596. {
  597.     register pset p;
  598.     register int i, base, *count;
  599.     register unsigned int val;
  600.     int weight;
  601.     pset last;
  602.  
  603.     count = ALLOC(int, A->sf_size);
  604.     for(i = A->sf_size - 1; i >= 0; i--) {
  605.     count[i] = 0;
  606.     }
  607.  
  608.     /* Loop for each set */
  609.     foreach_set(A, last, p) {
  610.     weight = 1024 / (set_ord(p) - 1);
  611.     for(i = LOOP(p); i > 0; ) {
  612.         for(val=p[i]&r[i], base= --i<<LOGBPI; val!=0; base++, val >>= 1) {
  613.         if (val & 1) {
  614.             count[base] += weight;
  615.         }
  616.         }
  617.     }
  618.     }
  619.     return count;
  620. }
  621.  
  622.  
  623. /*
  624.  *  sf_delc -- delete columns first ... last of A
  625.  */
  626. pset_family sf_delc(A, first, last)
  627. pset_family A;
  628. int first, last;
  629. {
  630.     return sf_delcol(A, first, last-first + 1);
  631. }
  632.  
  633.  
  634. /*
  635.  *  sf_addcol -- add columns to a set family; includes a quick check to see
  636.  *  if there is already enough room (and hence, can avoid copying)
  637.  */
  638. pset_family sf_addcol(A, firstcol, n)
  639. pset_family A;
  640. int firstcol, n;
  641. {
  642.     int maxsize;
  643.  
  644.     /* Check if adding columns at the end ... */
  645.     if (firstcol == A->sf_size) {
  646.     /* If so, check if there is already enough room */
  647.     maxsize = BPI * LOOPINIT(A->sf_size);
  648.     if ((A->sf_size + n) <= maxsize) {
  649.         A->sf_size += n;
  650.         return A;
  651.     }
  652.     }
  653.     return sf_delcol(A, firstcol, -n);
  654. }
  655.  
  656. /*
  657.  * sf_delcol -- add/delete columns to/from a set family
  658.  *
  659.  *  if n > 0 then n columns starting from firstcol are deleted if n < 0
  660.  *  then n blank columns are inserted starting at firstcol
  661.  *      (i.e., the first new column number is firstcol)
  662.  *
  663.  *  This is done by copying columns in the array which is a relatively
  664.  *  slow operation.
  665.  */
  666. pset_family sf_delcol(A, firstcol, n)
  667. pset_family A;
  668. register int firstcol, n;
  669. {
  670.     register pset p, last, pdest;
  671.     register int i;
  672.     pset_family B;
  673.  
  674.     B = sf_new(A->count, A->sf_size - n);
  675.     foreach_set(A, last, p) {
  676.     pdest = GETSET(B, B->count++);
  677.     INLINEset_clear(pdest, B->sf_size);
  678.     for(i = 0; i < firstcol; i++)
  679.         if (is_in_set(p, i))
  680.         set_insert(pdest, i);
  681.     for(i = n > 0 ? firstcol + n : firstcol; i < A->sf_size; i++)
  682.         if (is_in_set(p, i))
  683.         set_insert(pdest, i - n);
  684.     }
  685.     sf_free(A);
  686.     return B;
  687. }
  688.  
  689.  
  690. /*
  691.  *  sf_copy_col -- copy column "srccol" from "src" to column "dstcol" of "dst"
  692.  */
  693. pset_family sf_copy_col(dst, dstcol, src, srccol)
  694. pset_family dst, src;
  695. int dstcol, srccol;
  696. {
  697.     register pset last, p, pdest;
  698.     register int word_test, word_set;
  699.     unsigned int bit_set, bit_test;
  700.  
  701.     /* CHEAT! form these constants outside the loop */
  702.     word_test = WHICH_WORD(srccol);
  703.     bit_test = 1 << WHICH_BIT(srccol);
  704.     word_set = WHICH_WORD(dstcol);
  705.     bit_set = 1 << WHICH_BIT(dstcol);
  706.  
  707.     pdest = dst->data;
  708.     foreach_set(src, last, p) {
  709.     if ((p[word_test] & bit_test) != 0)
  710.         pdest[word_set] |= bit_set;
  711. /*
  712.  *  equivalent code for this is ...
  713.  *    if (is_in_set(p, srccol)) set_insert(pdest, destcol);
  714.  */
  715.     pdest += dst->wsize;
  716.     }
  717.     return dst;
  718. }
  719.  
  720.  
  721.  
  722. /*
  723.  *  sf_compress -- delete columns from a matrix
  724.  */
  725. pset_family sf_compress(A, c)
  726. pset_family A;            /* will be freed */
  727. register pset c;
  728. {
  729.     register pset p;
  730.     register int i, bcol;
  731.     pset_family B;
  732.  
  733.     /* create a clean set family for the result */
  734.     B = sf_new(A->count, set_ord(c));
  735.     for(i = 0; i < A->count; i++) {
  736.     p = GETSET(B, B->count++);
  737.     INLINEset_clear(p, B->sf_size);
  738.     }
  739.  
  740.     /* copy each column of A which has a 1 in c */
  741.     bcol = 0;
  742.     for(i = 0; i < A->sf_size; i++) {
  743.     if (is_in_set(c, i)) {
  744.         (void) sf_copy_col(B, bcol++, A, i);
  745.     }
  746.     }
  747.     sf_free(A);
  748.     return B;
  749. }
  750.  
  751.  
  752.  
  753. /*
  754.  *  sf_transpose -- transpose a bit matrix
  755.  *
  756.  *  There are trickier ways of doing this, but this works.
  757.  */
  758. pset_family sf_transpose(A)
  759. pset_family A;
  760. {
  761.     pset_family B;
  762.     register pset p;
  763.     register int i, j;
  764.  
  765.     B = sf_new(A->sf_size, A->count);
  766.     B->count = A->sf_size;
  767.     foreachi_set(B, i, p) {
  768.     INLINEset_clear(p, B->sf_size);
  769.     }
  770.     foreachi_set(A, i, p) {
  771.     for(j = 0; j < A->sf_size; j++) {
  772.         if (is_in_set(p, j)) {
  773.         set_insert(GETSET(B, j), i);
  774.         }
  775.     }
  776.     }
  777.     sf_free(A);
  778.     return B;
  779. }
  780.  
  781.  
  782. /*
  783.  *   sf_permute -- permute the columns of a set_family
  784.  *
  785.  *   permute is an array of integers containing column numbers of A which
  786.  *   are to be retained.
  787.  */
  788. pset_family sf_permute(A, permute, npermute)
  789. pset_family A;
  790. register int *permute, npermute;
  791. {
  792.     pset_family B;
  793.     register pset p, last, pdest;
  794.     register int j;
  795.  
  796.     B = sf_new(A->count, npermute);
  797.     B->count = A->count;
  798.     foreach_set(B, last, p)
  799.     INLINEset_clear(p, npermute);
  800.  
  801.     pdest = B->data;
  802.     foreach_set(A, last, p) {
  803.     for(j = 0; j < npermute; j++)
  804.         if (is_in_set(p, permute[j]))
  805.         set_insert(pdest, j);
  806.     pdest += B->wsize;
  807.     }
  808.     sf_free(A);
  809.     return B;
  810. }
  811.